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, Path Integral Quantum Monte Carlo simulation is used to study thermodynamic properties and 

' a phase diagram of 2D quantum Josephson array, described by 2+1 XY model. The helicity and 

t-H \ vorticity moduli, correlation function of phases and other characteristics of the system as functions 

» I . of quantum parameter q and temperature T are studied (q — 2e/V J Co, J is the Josephson coupling 

constant, Co is the intragrain capacitance). Quantum fluctuation induced superconductor - normal 
phase transition is studied in detail through the use of behavior of above-mentioned quantities. No 
discontinuous or reentrant phase transition in q — T plane is found. Analysis of the vorticity and 
l/") ' the renormalized coupling constant leads to the conclusion that the whole line of phase transition 

is of Kosterlitz - Thouless type. 

Q . I. INTRODUCTION 

; 

!_i ■ Ordering phenomena in 2D Josephson arrays has elicited considerable theoretical and experimental interest (see jl[] 
Oh| - |?J and references therein). The phase transition in classical system when capacitances are neglected is generally 
believed to be of Kosterlitz - Thouless type. Taking into account intra- and intergrain capacitances leads to a 
variety of interesting phenomena, which are governed by the competition between the Josephson tunneling and the 
charging energy of the grain. This charging effect leads to the breaking down (at appropriate sizes of granules) of 
superconductivity even at zero temperature. 

1 I A major outcome of the simulation of periodic array of ultrasmall granules |^| was the observation of a transition 

to a new coherent state. At sufficiently small capacitances it was shown to be two " ordering" phenomena: the system 
undergoes at higher temperatures Kosterlitz - Thouless transition from the disordered phase into the superconducting 
one, and upon further decreasing of temperature a quantum induced first order transition to another superconducting 
phase with finite but diminished superfluid density. 

It is important to analyze a possible modification of mechanism and type of the phase transition when quantum 
effects begin to play an important role. Because the critical behavior of the quantum system is similar to that 
of the thick classical anizotropic film (and is described in terms of 2+1 XY model), the Kosterlitz-Thouless phase 
transition scenario may be in principle replaced (for the finite size system) by behavior, inherent to 3D anizotropic 
(f-) I superconductor ||, 0M. By introducing a Hikami-Tsuneto (HT) length scale defined by the ratio of the inter (Jj_) 
(T) . and intra(J||) - plane coupling constants |n]] Pht = y/J±J J\\ = P/ilPJ) (P - is the thickness of the anizotropic 
classical film, (3 = l/k^T) it can be easily seen, that at distances \r*i — fj\ < P/Pht the behavior of the system has 
3D character with a corresponding correlation function Q and vortex characteristics, while at \ fi — fj\ > P/ Pht the 
system has 2D behavior. Thus, it is interesting to analyze topological excitations responsible for phase transition in 
the region of strong quantum fluctuations. 

The purpose of this Letter is ab initio studying properties and the phase diagram of the 2D quantum Josephson 
array. By path integral Monte Carlo simulations (PIMC) we will calculate different physical properties (helicity 
and vorticity moduli, correlation function, and vortex moment) on the plane {q,T}, where the quantum parameter 
q = 2e/v 'JCq is responsible for the strength of zero - point fluctuation of phases, T = kbT/J is the dimensionless 
temperature of the system. From the analysis of above-mentioned quantities we will show, that at distances larger 
than Pj Pht the system experiencing 2D Kosterlitz-Thouless - like phase transition with vortex unbinding and thus 
can be regarded as 2D classical system with a coupling constant renormalized by quantum fluctuations. 

X 

^ . II. MODEL AND CALCULATING PROPERTIES. 

At sufficiently low temperatures fluctuations of moduli of the superconducting order parameter can be neglected 
and the Hamiltonian of the 2D system of TV 2 granules can be written as 

2e 2 

H = — V(rii - n ) 2 + J V (1 - cos (<& - ipj)) 

On ^""^ . 7 

i <t,]> 

The first term of the Hamiltonian is related to the Coulomb charging energy (Co is the intragrain capacitance; We 
suppose that intergrain capacitances C\ 3> Co and their effect may be neglected) . fii — uq stands for the deviation of 
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the number of effective bosons (Cooper pairs) from its equilibrium value no =< hi >. Josephson coupling energy 
J in the second term is assumed the same for all nearest - neighbor pairs < i,j > of islands. Periodical boundary 
conditions are used. Phases ifi of the order parameter are chosen to be cyclic variables: (fi £ [0, 2ir), i — 1 . . . TV 2 . 

If the relative fluctuations of the number of Cooper pairs are small then the particle number operator hi — hq can 
be chosen as one conjugate to the 'phase' operator |ll],[l^]: hi — uq = jd/d<fi. This lead to the Hamiltonian of 2D 
quantum XY model: 
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After the simple discretisation procedure for path integral, in the so - called primitive form of PIMC (see e.g. |JL3]), 
one can represent properties of an initial 2D quantum system by fictitious 3D classical system formed by the P - 
times multiplication of the initial one. It is described by Boltzmann partition function w with an effective potential 
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E — - — h-2 + tp E c 1 - cos (^ - $)) 
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V? = • (2) 

an integer n* provides the minimum of |</?f +1 — <f% — 2im*\ at some fixed , and approximately takes into 
account that phases tp? are cyclic variables. An equilibrium value of some observable O then can be found as: 

/ . . . 7o(<pl<4, . . . pSUM^S ■ ■ • ■ ■ ■ d^ 1 

O) - (0) p - ° ° 2 , 2 . (3) 



In principle we can substitute the fictitious classical system in hand into the anizotropic 3D XY one with the coupling 
constants Jii = J/P, J± = P/{q 2 J(3 2 ). The anizotropy of this system is defined by the ratio J±/J\\ — 2/n, were 
the parameter t\ <C 1 is responsible for the accuracy of the primitive discretisation procedure (see below). So, one 
can imagine, that the behavior of the anizotropic 3D XY system N x N x P with a constant anizotropy is studied at 
different points {q, T} of the phase diagram and, consequently, with different thickness P. 

Properties of the system were explored in approximately 150 points belonging to the region 0.1 < q < 2.6, 0.02 < 
T < 1.4. We used the multigrid algorithm (unigrid modification Q was used). More precisely, we used the Ww cycle 
(i.e. the W cycle with 1 presweep and no postsweeps) in the "classical region" q < 0.5 while the Vio cycle was found 
to be more efficient at q > 0.7. The Trotter number P at each point of phase plane was fitted in order to the relative 
discretisation error of Feynman integral being smaller that 3%. This type of error is determined by the discretisation 
parameter t\ = 2q 2 /T 2 P 2 , its value was 0.001 < T\ < 0.07 at all considered points. Note, that the errors related with 
partial neglecting of the periodical boundary conditions for temperature density matrix (the appropriate parameter is 
T2 = 2 cxp (— 2n 2 T P / q 2 ^j) are small with such Trotter numbers and are of order of 2% of the calculated mean values. 

The following quantities where measured at all above-mentioned points of the {q, T} plane: 

1) The modified Lindeman ratio. 

This quantity is the generalization of the well-known Lindeman's rms displacement of particles (phases). The 
Lindeman criterion states that solid melts when the rms displacements become larger than some universal part 5 C 
of the interparticle distance (5 C is generally assumed to be » 0.15). For 2D system this value is known to diverge 
logarithmically with the system size. This does not enable us to take advantage of this universality in determining 
the temperature of phase transition. Instead of rms displacement one can use the modified Lindeman ratio Jl5| ]: 

Sl = \2tP E LVi - <Pi\\-*,*) ) ( 4 ) 

\ <i,j> > l I 

The sum in the formula (Q) is taken over such pairs of granules < i, j > that the vector drawn from the granule i to 
j is equal to (0, 1) or (1, 0). In view of periodic boundary conditions the number of such pairs is equal to 2iV 2 . Here 
and henceforth [f\ j a b ^ denotes the reduction of / to the interval [a, b). 
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2) The helicity modulus. 

This quantity gives the information about the "rigidity" of the system perturbed by imposed gradient 5<j){m,n) = 
kf(m, n) — k x m + k y n , r(m, n) £ [1, AT] x [1, N]. 

Helicity modulus Jl6| (see also ||T^| ) can be defined as "f a b = -jjjv gf g kb where F — —kbTln(Z) is free energy 

of the system. It is worth-while to note, that the classical expression (see e.g. |lq| ) for 7 is not correct for quantum 
system. The reason is the same as that pointed out in the review [ p"3| when calculating superfluid density: operators 
dH/dk and e~@ H do not commute. After some algebra, the following expression for the helicity modulus of the 
quantum XY model can be obtained JsJ : 

lab = jSab, 
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Jp ( cos(tp°(m + l,n)-ip°(m,n))\ - 

\m,n=l / p 

sin {<p°{m + 1, n) - cp°(m, n)) sin {<p p (k + 1,1)- ip p (k, I)) \ (5) 
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On assuming that ip p (m, n) — (p°(m, n), one can obtain the well-known classical expression |l 
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7ds = -TTn { cos (cp (m + 1, n) - tp (m, n)) 
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sin (ip°(m + 1, n) — ip°(m, n)) 
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Being calculated, the "classical helicity modulus" can be regarded as a low bound to the true one (|5|). It occurs that 
the classical expression gives rather a good approximation of (|^) almost for all the superconducting region of the 
{q, T} plane excluding q > 2 and small temperatures T < 0.05 where the number of effectively uncorrclatcd planes 
P/Pht is large enough and the imaginary - time correlator in (||) is far from constant. 
3) The space correlation function. 

The space correlation function g(l) of phases can be introduced as: 



5(0 = ( — ^ cos ('Pi ~ <P°) ) ( 7 ) 




where n\ is the number of granules at the distance of I (0 < I < [yj — 1) from some fixed one <p c . 

4) The density of open vortices. 

If there are 2n open vortex threads in the system [0, TV] x [0, N] x [0, P] (see e.g. then the density of vortices 
p v is defined as: 

* - W (8) 

5) "Vorticity modulus". 

This quantity is the measure of the response of the system to an isolated vortex (open vortex thread in our case) 
introduced into the system and may be defined as p3] : 

Fdap - F P 

VDAP - P = In TV (9) 

Here Fp is free energy of the system under periodic ('P') boundary conditions when the net topological charge is 
equal to zero and there always appear only equal numbers of vortices and antivortices. Fdap is free energy of the 
system under diagonally - antiperiodic ('DAP') boundary conditions (see below) when there is at least one excess 
vortex in the system. 

Variational upper and lower bounds on free energy differences (bl) can be obtained via Bogoliubov - Gibbs variational 



method (see e.g. 21 1). To find upper bound, for example, the coupling constant at the boundary J t & , and consequently 
the Hamiltonian, is changed slowly over the course of the simulation (where the parameter t = . . . 2M marks the 
Hamiltonian updates) from its initial value Jg = J to zero — (free boundary conditions), system being at 
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periodic boundary conditions. At this the very moment t — M the type of boundary conditions is switched to DAP 
i.e. 

^(m, 0) = [tt + y»{N - m +l,N)\ [0>!hr) , <f?(m, N + 1) = [tt + ip p (N - m + 1, 1)J [0 , 27r) , 

^(0, n ) = [tt + tp p (N, N - n + 1)J [0>27r) , ^(A + 1, n) = + (p p (l,N -n+ 1)J [0i27r) . 

Subsequent increasing of the coupling constant on the boundary up to its initial value J^ M — J led system to the 
final state t = 2M with at least one free vortex presented. 

At each step t the system was heated for ~ 400A 2 Monte Carlo (MC) steps (one MC step is the V cycle or 
the W cycle). Measuring of the work 5A t of altering the coupling constant from J(t) b to J(t + l) b required some 
more ~ 400A 2 MC steps. The energy difference caused by alteration of the coupling constant at the boundary 
J(t) b — > J(t + l) b contributes to the work as: 



SA t = (H(J b +1 ) - H(J b )) Jb = ( J b +1 ]T(1 - cosfe - - J b ^(1 - cosfo - <pj)) 



The upper estimation of free energy difference (Q) is calculated as the total work required for changing of boundary 
conditions from P to DAP: 



F DA p -F P <Y,SA 



2M 

P^DAP 
t 

t=0 



The calculation of the lower bound of free energy differences required the analogous transformation from the DAP 
to P. Similarly, we obtain: 



F D ap-Fp > -J2 6A i 



2M 

DAP^P 



t=0 

In present work, we found that a good choice of slowly varying coupling constant at the boundary is 

J b = J cos 2 (nt/2M). 

The concrete number of Monte Carlo steps for heat and measure as well as the total number 2M of Hamiltonian 
updates where fitted to provide the accuracy (it can be estimated as the difference between the upper and the lower 
bounds of the vorticity modulus) needed for the estimation @. 

There is another way of looking at the vorticity modulus. Let us assume that the system is placed to a magnetic 
field H such that: 

, . . a % &os , ^ Inch 

(A x , A y ,A z ) = Mxa+y2) {-V, x, 0), $ = — 

(H X ,H V ,H Z )= 8 ^ + ^ (0,0,1) 

TT 

Let the origin of the frame of reference be in the center of the array. Consider the free energy of the system as a 
function of the strength s of the field. In the limit of small s one can write 

F(s)=F(0) + ^ S 2 + O( S 3 ) 

The linear response JV = d 2 F(s)/ds 2 \ s _ of the system to an isolated vortex at the origin takes the form: 

/ N 
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V=( ^2 {A 2 x (m,n) cos (tp a {m + l,n) - ip°(m,n))+ A 2 (m, n) cos (<p°(m, n + 1) - <p°(m, n))}) p - 

\m,n— 1 
P-l N 

^ {A x (m, n) sin (<p p (m + 1, n) — <p p (m, n))+ A y (m, n) sin (^ p (m, n + 1) — (p p (m, n))}] 2 ^ , (10) 



p—0 m,n=l 
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A x (m,n) = arctan (z m ,„/y m ,„) - arctan {x m+1 . n /y m . n ), 
A y (m,n) = arctan (y m ,n+i/x m , n ) ~ arctan (y m , n /x m ,n), 

where (x m ^ n , y m ,n) stands for the position of the granule (m,n) in the particular frame of reference. 

Analogous expression was used in studying the phase diagram of the Heisenberg antiferromagnet [p2| . As have been 
pointed out, the V contains both a core energy E c and a part which is proportional to In AT: 

V(N) = E c /J + v H \nN 

The size dependent part can be extracted then by using the results obtained for systems of size N\ and N 2 [f22j : 

V(N 2 )~V(N 1 ) 
VH = ln(N 2 / Nl ) 



III. RESULTS AND DISCUSSION. 

Figures 2a,b show behaviour of the modified Lindeman ratio (0) in the explored region of parameters {5, T}. From 
the results of calculations one can see that neither in the 'classical' region (q < 1) nor in the quantum one, this quantity 
does not have any points of abrupt changes that could have revealed the line T c (q) of phase transitions. But it is hoped 
that the values Sf at the point {q c , T c } of phase transition may be universal for all 2D system of the same symmetry. 
For example, in the case of classical system (q — 0) we have = <5i (0, 0.89) sw 0.81, £3 = £3(0, 0.89) sw 1.06. It should 
be particularly emphasized that the abovementioned quantities Si have another values 5[ at the point {q, T} = {q c , 0} 
of quantum phase transition (this universal values, if any, are the same for analogous quantum system). 

Disappearance of the helicity modulus testifies to the disordering of phases in the system. Let us analyze the 
behavior of 7(9) (see Fig. 3a). As the quantum parameter q increases, the value of 7 decreasing slowly up to some 
critical q* rapidly falls to zero in a narrow interval (q* , q* + Aq) . Analogous behavior of this quantity takes place 
along all lines q = const. Really, at q — there is a sharp drop of the quantity j{T) at T* w 0.9, associated with the 
Kosterlitz - Thouless phase transition in the classical system. At Figure 3b examples of such classic - like behavior of 
7(T) for different values of the quantum parameter q are presented. From the Figures 3 one can see, that the curves 
7(g) and j(T) are very gently sloping in the region q > 1.0. 

In the Fig. 3a (solid diamonds) the results of calculations of the helicity modulus at q — 0.5 and at low temperatures 
T < 0.1 are given. In this the very region of temperatures, at T < 0.03 authors of the work j|] marked a discontinuity 
in helicity modulus and specific heat as functions of temperature. This phenomena was associated with the first order 
phase transition due to quantum fluctuations of phases. Presented results correspond to the system N x N = 10 x 10 
(as well as in the ||), but our discretisation errors are much smaller: Trotter number P was P|t=o.02 = 128 (n ~ 0.07 
compared to t\ ~ 0.25 in the ||). We see that within the limits of statistical errors any disordering phenomena are 
absent in this region {q,T}. No reentrance (or discontinuity) are observed (see Fig. 3b) also at larger values of the 
quantum parameter q, i.e even when a strength of quantum fluctuations is rather great. It is worth-while to note, 
that the classical expression for the helicity modulus (^|) is incorrect at this region of dimensionless parameters of the 
system. The deviation of the value of 7 calculated by the adequate expression (|B|) from that calculated by classical 
expression may have the order of 7 both at small temperatures T < 0.05 (at q ~ 0.5) and at large values q > 2 of the 
quantum parameter, to lead to the appreciable altering of the phase diagram. This suggestion is supported by results 
presented at Figures 3a,b. 

As pointed out in || , the helicity modulus (the superfluid density) tends to a finite value of order one as T — ► 
(we have 7| r _ ~ 0.87 at q ~ 0.5). So, we can apply the spin wave consideration to the anizotropic system 
N x N x P in this region of temperature. Through the appropriate scale transformation we obtain the izotropic 
system N x N x P/Pht with coupling constants J11 = J± = Pht/P, Pht = PT/q. Note, that the thickness of this 
system (T = 0.02, q = 0.5, N = 10) is three times larger than its width: P/Pht = q/T ~ 25 > N and hence we are 
to observe 3D - ID crossover behavior at this parameters of the system. 

The quantum N x N Josephson array is known to have 3D - like behavior below the line T — q/N ||, In the 

above region we are to observe 2D behavior and, hence, the boundary of the global superconductive state in q — T 
phase diagram (T > q/N) is expected to be a line of Kosterlitz - Thouless transitions with temperatures T c — T c (q) 
defined via. the renormalized by quantum fluctuations coupling constant J{q). 

Insufficient accuracy of calculations and small system sizes does not permit, as a rule, to check for agreement 
between theoretical Kosterlitz - Thouless critical exponents and calculated ones. This implies that the conclusion 
about a type and a character of the phase transition may be done on the basis of some other specific properties. 
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The vorticity modulus is one of them for topological Kosterlitz - Thouless phase transition and is proportional to 
excess free energy due to an excess vortex. In terms of the 2+1 XY model we may consider that this quantity is the 
measure of the ability of the N x N x P system to form an isolated unclosed vortex thread. From the definition of 
this quantity (^) it follows that the vorticity modulus vjjap-p being positive in the ordered phase, should vanish at 
the disordered one. Dependencies VDAP-p(T)\ q=const for the system N x N = 12 x 12 and different values of the 
quantum parameter are presented in Fig. 4a (the dependencies vjj (T) and vu (q) , see ([ll]) , are given at Figures 5a, b) . 
Shapes of curves VDAP-p(T)\ q=const for classical (q — 0) and quantum regions are very much alike. The tendency of 
decreasing the phase transition temperature with increasing q correlates with the behavior of the vorticity modulus: 
the temperature T v at which the vorticity changes sign, i.e. it is advantageous to born an excess vortex thread, is 
shifted to the lower temperatures as quantum parameter q is increased. This observation is in agreement with the 
results of calculations of the helicity modulus vs. temperature at different q (see Fig. 3b). The dependence of the 
vorticity vs. q at T — 0.5 is given at Fig. 4b. As was noted by authors of |po| , the temperature T v decreases only 
weakly with the system size. This is the reason, by which we may use data on calculating vorticity modulus for the 
system 12 x 12 together with other results for the system 20 x 20. 

The point {q v ,T v } at which the vorticity modulus (|) changes sign may be considered as one at which the mean 
separation of free vortices matches the system size So it would be interesting to compare points {q v ,T v } with 
those of great variance of the correlation length of phases. Let us suppose that an empirical formula for the correlation 
function (Q) g(r) ~ exp (— r/£) takes place. Experimental data on £ as a function of q (or T) is presented at fig. 4b 
(Fig. 4a). It may be seen, that points {q v ,T v } are belongs to the regions of considerable increase of the empirical 
parameter £ which may be associated then with the correlation length of the system. 

An analysis of the helicity and the voricity moduli make it possible to consider the point {q,T} = {2.2,0.5} = 
{<7 C ,T C } as the point of phase transition. Let us explore the region q ~ 2, T ~ 2 of strong quantum fluctuations 
near this point. The corresponding 3D effective izotropic system in this region has dimensions N x N x P/Pht ~ 
20 x 20 x 4 and can be regarded as 2D one at distances |r, > 4. Particles of this 2D system, interacting via the 
coupling constant J(q, T) renormalized by quantum fluctuations, can be associated with open threads (of appropriate 
topological charge). To determine the renormalized coupling constant J = J(q,T) of the 2D system of charges at 
some point {q, T} one can analyze the density p v (||) of open threads in the system. When p v << 1 one can write 



P 



-(l- N2-2^l3,J(q.T)^ 



After p v (q c ,T c ) was being calculated, we compared the experimental renormalized coupling constant J(q c ,T c ) with 
that, obtained by Kosterlitz - Thouless criteria: J c = T c /0.93 » 0.54. Rather a good conformity between both the 
results testify to the correctness of the Kosterlitz - Thouless picture of phase transition in this quantum region (more 
detailed analysis of quantities described the vortex structure of the system will be given elsewhere) . 

Another argument is the universal jump of the helicity modulus in this the very point {q c ,T c }. The point qj of the 
universal jump is that of intersecting of the line 2T/7r| T=0 5 — 1/tt with the plot of the helicity modulus "f(q)\ T=0 5 
(see Fig. 3a) and is equal to qj — 2.15 ± 0.05. 

The resulting phase diagram (see Fig. 1) is in sufficiently good agreement with that calculated from Kosterlitz - 
Thouless scenario through the use of self-consistent harmonic approximation (SCHA) p],p4|]. 

To conclude, we have analyzed properties of 2D Josephson array at different quantum parameters q and temperatures 
T. The phase transition in all the region (above the region of crossover) have been found to be of Kosterlitz - Thouless 
type. No reentrance or discontinuity phenomena have been found. 
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Fig. 1. 

The phase diagram of 2D quantum Josephson array. (1) Results of present calculations. (2) Results of Ref. (3) 
SCHA renormalized Kosterlitz - Thouless temperatures p3 , 

Fig. 2a. 

Modified Lindeman ratio Si as a function of temperature T. Data are connected to guide the eyes. If not presented, 
the error bars are within the size of the data point. 

Fig. 2b. 

Modified Lindeman ratio Si as a function of quantum parameter q at T = 0.6. 
Fig. 3a. 

The dependence of the helicity modulus 7 (and 7 C ; S ) vs. quantum parameter q at temperature T = 0.5. N x N = 
20 x 20. The dependence — = 1/n (see in text) is shown with the help of a dashed line. 

Fig. 3b. 

The dependence of the helicity modulus 7 (and 7 C ; S ) vs. temperature T at different quantum parameters q: Insert: 
3 = 0.5(^x^ = 10x10); 

Fig. 4a. 

The vorticity modulus vdap-p @ as a function of temperature T at different values of the quantum parameter. 
N x N = 12 x 12). 

Open squares are the empirical coefficient £ = £(T) in the expression for the space correlation function of phases 
g(r) ~ exp(-r/£); g = 0.7, N x N = 20 x 20. 

Fig. 4b. 

The vorticity modulus vuap-p vs. quantum parameter q at T = 0.5, N x N = 12 x 12. Open squares are the 
empirical coefficient £ = £(g) (see in text) in the expression for the space correlation function of phases (^) at T = 0.5, 
N xiV = 20 x 20. 

Fig. 5a. 

The vorticity modulus vh ( pT| ) as a function of temperature T at q = 1.5. 
Fig. 5b. 

The vorticity modulus vh ( |lT| ) as a function of q at T = 0.6. 
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